Create algorithm for finding direction of fibers

This algorithm uses skewness in fourier transform to calculate direction of fibers in image. It does so by fitting a line to points above a given threshold in the fourier spectrum.


In [1]:
%pylab inline
from skimage.io import imread
from skimage.external import tifffile #only in v0.11


Populating the interactive namespace from numpy and matplotlib

Turn of matplotlib interpolation, set gray colormap and figure size


In [2]:
plt.rcParams['image.interpolation'] = 'none'
plt.rcParams['image.cmap'] = 'gray'
figsize(4,4)

Open a console to do temporary work


In [186]:
%qtconsole

Create a sample image


In [3]:
size = 256
img = np.zeros((size,size), dtype=np.uint8)
t = linspace(start=0, stop=50*pi, endpoint=False, num=size)
x,y = meshgrid(t, t)
img[:,:] = 127 + 127*sin(x+0.2*y)
imshow(img);


Do fourier transform and find a suitable threshold


In [4]:
F = log(abs(fftshift(fft2(img))))
thresholded_F = np.copy(F)
threshold = F.max() * 0.25
thresholded_F[F<threshold] = 0
f, axs = subplots(ncols=2, figsize=(8,4))
axs[0].imshow(F)
axs[1].imshow(thresholded_F);


Find threshold based on number of pixels (highest 1‰ of them)


In [5]:
from skimage.exposure import cumulative_distribution
cdf = cumulative_distribution(F)
plot(cdf[1], cdf[0])
limit = np.where(cdf[0] > 0.999)[0].min()
threshold = cdf[1][limit]
plot([threshold, threshold], [0,1]) # vertical line
threshold


Out[5]:
7.0729400447769724

In [6]:
thresholded_F = np.copy(F)
thresholded_F[F<threshold] = 0
f, axs = subplots(ncols=2, figsize=(8,4))
axs[0].imshow(F)
axs[1].imshow(thresholded_F);


Find a line which fits the data


In [7]:
# points
xx,yy = np.where(thresholded_F)
# y = mx + c -> matrix equation: y = Ap where A = [[x 1]] 
# and p = [[m], [c]]
A = np.vstack([xx, np.ones_like(xx)]).T
p = np.linalg.lstsq(A, yy)[0]
# polynomial generator
poly = np.poly1d(p)

# plot it
line = [poly(0), poly(F.shape[1])]
plot([0,F.shape[1]], line)
xlim(0, F.shape[1])
ylim(F.shape[0], 0);



In [11]:
# assume the line is passing through center
angle = arctan(p[0]) / pi * 180

print(angle)

arctan(-0.2/1) / pi * 180


6.87033788935
Out[11]:
-11.309932474020215
  • $\theta = 0$ should be a horizontal fiber
  • positive y axis is down ==> angle goes clockwise

==> we should have angle $\approx$ +80$^\circ$

  • polyfit is turned 90$^\circ$ right in regard to the fiber

Lets make degrees go counterclockwise instead (it's hard to "think around" convention):

$\theta = tan^{-1}(\frac{-y}{x})$

And convert to 0-180 range:

$\theta^` = \theta \% 180$


In [12]:
# turn left 90 degrees
angle = arctan(-p[0]) / pi * 180 + 90
# range 0-180
angle = angle % 180
angle


Out[12]:
83.129662110646279

Make a function of it all


In [13]:
def ft_line_fit_angle(img, threshold=0.999, debug=False):
    """
    Calculate preferred orientation in image with a line fit in FT.
    
    Parameters
    ----------
    threshold : float
        Percentage of pixels to include in FT for calculating 
        threshold. 0.999 * 512**2 = 262 pixels
        
    Returns
    -------
    float
        Angle
    """
    # FT power spectrum
    F = np.abs(fftshift(fft2(img)))
    # do not calculate log(0)
    F[F!=0], F[F==0] = log(F[F!=0]), log(F[F!=0].min())
    # threshold
    cdf = cumulative_distribution(F)
    limit = np.where(cdf[0] > threshold)[0].min()
    threshold_value = cdf[1][limit]
    F = F > threshold_value
    # points
    y,x = np.where(F)
    # special case (trivial solution)
    if x.min() == x.max():
        # we have a vertical line
        angle = 90
        b = [0, 1]
    else:
        # y = mx + c -> matrix equation: Ab = y where A = [[x 1]] 
        # and b = [[m], [c]]
        A = np.vstack([x, np.ones_like(x)]).T
        b = np.linalg.lstsq(A, y)[0]
        # calculate angle (assume line goes through center)
        angle = (90 - arctan(b[0]) / pi * 180) % 180

    # show image, FT and fit
    if debug:
        f, ax = subplots(ncols=2, figsize=(8,4))
        ax[0].imshow(img)
        ax[1].imshow(F)
        # add calculated line
        # polynomial generator
        p = np.poly1d(b)
        height, width = img.shape
        if angle != 90:
            line = ([0, width], [p(0), p(width)])
        else:
            line = ([width//2, width//2], [0,height])
        ax[1].plot(*line)
        ax[1].set_title('computed angle: {:3.0f}'.format(angle))
        ax[1].set_xlim(0,width)
        ax[1].set_ylim(height,0)

    return angle

Test it with another orientation of "fibers"


In [14]:
angle = 4
dx = sin(angle/180*pi)
dy = cos(angle/180*pi)
img = 127 + 127*sin(y*dy+x*dx)
ft_line_fit_angle(img, debug=True)


Out[14]:
90

A range of angles


In [484]:
for angle in range(0,180, 10):
    dx = sin(angle/180*pi)
    dy = cos(angle/180*pi)
    img = 127 + 127*sin(y*dy+x*dx)
    ft_line_fit_angle(img, debug=True)


Lets try with an image


In [539]:
from itertools import product
img = imread('mp.png')
bs = 100 # block size
yy,xx = img.shape
blocks_y, blocks_x = yy//block_size, xx//block_size
count = 0
for j,i in product(range(blocks_y), range(blocks_x)):
    x = j*bs
    y = i*bs
    temp_img = img[y:y+bs, x:x+bs]
    if temp_img.shape[0] < 50 or temp_img.shape[1] < 50:
        continue
    ot = threshold_otsu(temp_img)
    if (ot < 4 # threshold above noise
        or (temp_img > ot).sum() < 20): # at least 20 pixels
        continue
    count += 1
    if count > 10:
        break

    #figure(figsize=(2,2))
    #imshow(temp_img)
    #title('{:2.2f} {:2.2f}'.format(temp_img.mean(), threshold_otsu(temp_img)))
    ft_line_fit_angle(temp_img, threshold=0.9, debug=True)


Iterate over blocks of the image, and make a histogram of it all


In [14]:
from itertools import product
from skimage.filters import threshold_otsu
img = imread('mp.png')
bs = 100 # block size
blocks_y = img.shape[0]//bs
blocks_x = img.shape[1]//bs

# spread spare pixels
bsy = img.shape[0] // blocks_y
bsx = img.shape[1] // blocks_x

h = np.zeros(180) # histogram
for j,i in product(range(blocks_y), range(blocks_x)):
    x = j*bsx
    y = i*bsy
    temp_img = img[y:y+bsy, x:x+bsx]
    
    if temp_img.shape[0] < 50 or temp_img.shape[1] < 50:
        # small image
        continue
    
    ot = threshold_otsu(temp_img)
    if (ot < 4): # threshold below noise
        #or (temp_img > ot).sum() < 20): # at least 20 pixels
        continue

    angle = ft_line_fit_angle(temp_img, threshold=0.9)
    angle = round(angle)
    h[angle] += 1

Plot histogram


In [15]:
figsize(8,4)
plot(h);


Basic algorithm seems ok, but we see that angles are distributed around 90$^\circ$, which seems odd. Find cause and optimize the algorithm in optimize ft line fit.ipynb.